Topology and symmetry breaking in ABC trilayer graphene 
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The effects of topology and electron-electron interactions on the phase diagram of ABC stacked 
trilayer graphene (TLG) at the neutrality point are investigated within a weak coupling renormal- 
ization group approach. We find that the leading instability of TLG with only a forward scattering 
density-density interaction is towards a mirror-breaking gapless state. Addition of a small, but finite 
back scattering favors gapped phases, allowing us to make connections to the existing experiments 
on TLG. We identify a fundamental symmetry difference between TLG and bilayer graphene (BLG) 
which is responsible for disfavoring nematic states in TLG under the same conditions that favor 
nematic states in BLG. 



Interest in the theory of metal-insulator transitions in 
divalent materials [l[ has undergone a resurgence as of 
late in no small part due to the discovery of graphene 0- 
15]. There is growing experimental evidence for electron- 
electron (e-e) interaction-induced excitonic instabilities 
in AB stacked bilayer (BLG) I5l4ll| and the ABC stacked 
trilayer graphene (TLG) [121 ]. Theoretically, based on 
symmetry alone TLG is qualitatively different from 
single layer graphene (SLG) and BLG. The valence band 
(VB) and the conduction band (CB) in SLG and BLG are 
members of the same irreducible representation (IR) of 
the symmetry group at the ±K = ±-^=^x points in the 
Brillouin zone. On the other hand, the VB and the CB 
of TLG belong to different IR's such that the interband 
matrix clement of the momentum operator vanishes at 
±K. 

SLG, BLG, and TLG are further distinguished in the 
topological character of their VB and the CB Bloch states 
as they wind along a closed k-space loop around the ±K 
points. Their winding numbers, which are 1, 2 and 3 re- 
spectively [l3, 14 1, can be related to the experimentally 
measured steps in the quantized Hall conductivity as the 
carrier concentration is tuned from hole- to electron-like 
SEE 13- This non-trivial topology guarantees at least 
three Dirac Fermions (DF) near, but not exactly at, each 
of ±K of BLG and TLG, placed along the 3-fold symme- 
try lines in the BZ. Since one can adjust the microscopic 
parameters of the model without lowering the symmetry 
so as to either rotate the three DF's or to introduce an 
additional six DF's near each ±K, the Dirac degenera- 
cies in TLG cannot be established based on symmetry 
alone [17j ]. This raises the possibility of zero tempera- 
ture phase transitions between different, but symmetry 
indistinguishable, thermodynamic phases. 

In this work we systematically investigate the inter- 
play of such topological effects and the e-e interactions in 
TLG. We employ a powerful theoretical machinery based 
on symmetry and renormalization group (RG), which 
leads us to the basic structure of the general phase dia- 
gram shown in Fig. [T] together with the list of, and the 
conditions for, the most dominant symmetry breaking 



(SB) phases. As a biproduct, our symmetry analysis al- 
lows us to readily explain why 3-fold rotational SB gap- 
less states in TLG considered recently within Hartree- 
Fock approximation are indeed unfavorable 1181 . unlike 



in BLG msi, 

while noting that there are 3-fold rota- 
tionally symmetric gapless states, not considered in Ref. 
Hil which break the mirror symmetry of the lattice and 
are stable in the limit of pure forward scattering. Never- 
theless, small, but finite backscattering (i.e., intcrvalley 
scattering) makes gapless states unstable and promotes 
gapped phases. Our theory focuses from the outset on 
the most important low energy degrees of freedom and 
extends si gnif icantly beyond the recent functional RG 
treatment [23[ that did not take into account RG rel- 
evant symmetry allowed terms and treated only up to 
nearest neighbor lattice interactions. 



Because the experimentally observed 2J] energy (few 
meV) scales associated with any ordering are small com- 
pared to the relevant band splitting scales (few 100's 
meV), the theory of the TLG needs to take into account 
only the modes which lie close to the Fermi level. The 
effective low energy Hamiltonian written in the vicinity 
of the ±K points is 



H = Hn 



Hi, 



(1) 



Setting h = kg = 1, the symmetry allowed terms in the 
non-interacting Hamiltonian are 



1**, (2) 



where V = + d^^) 1 ^ + ( A + ^k 2 )!^, 

h^' = R x k x {k 2 x - 3kl)(R a T 3 t 2 + t 3 oi) and = 
i?2^ y (3fc 2 — fc,y)l2(T2. The Pauli matrices r and a 
act in the valley and layer space respectively. For 
each spin projection, the Fermi annihilation operators 
have been put into a four component vector ipk a = 
( C k+K CT > c^ +K(T , c£_ K(T , cf_ K(T ) T and |k| < A where A < 
|K|. A and B refer to two undimerized sites on layers 
1 and 3 (see Fig. [3] inset). There is some uncertainty in 



2 



the actual values and signs of the coefficients m* and A 
[H-[27|. In this Letter we choose m* > and sweep over 
a range of both positive and negative A's. A situation 
where m* < can be mapped onto the m* > problem 
by simultaneous change in sign of m* and A. Therefore, 
the results presented here cover every possible arrange- 
ment of parameter signs. 

For A > A c i = there are three anisotropic DF's 
along three symmetry lines, and we label this state 3 + 
(see Fig. [1]). For A < A C 2 < A c i there are also three 
DF's but rotated by 180° relative to 3+. We label this 
state 3~. For A C 2 < A < A c i, and in the presence of 
particle hole symmetry, which we assume to hold from 
now on, Rq = 0, and 1/Sm* = 0, and the spectrum 
contains nine DFs; we label this state 9~ (see Fig.[T]) [28j |. 
Because at T = it is impossible to transition between 
any two of 3 + , 3~, and 9~ states without encountering 
a non-analyticity in the ground state energy, even in this 
non-interacting case, there are two T = quantum phase 
transitions. The transition at A c i is 2 nd order and the 
one at A c2 = -l/(54(i?i) 2 m* 3 ) is 3 rd order 

We find that, for arbitrarily small e-e interactions, the 
2 nd order phase transition, initially at A cl , is avoided 
and replaced by two continuous phase transitions into a 
spontaneously broken symmetry phase (see Fig. [1}. The 
nature of the SB phase depends on the type and the 
strength of the e-e interactions and we discuss it in de- 
tail later in the text. On the other hand, the 3 rd order 
phase boundary point at A C 2 turns into a critical line 
A C 2(s) for a finite range of weak interactions g, before it 
is terminated in quantum tricritical points. 

We construct Hint using the irreducible representa- 
tions (IR's) of the space symmetry group P3ml [30(. The 
most relevant scattering processes are described by the 
quartic contact terms. These are the products of two bi- 

lincars gMN (Vi(r)M Q ^/3(r)) (vi(r)^ApVv(r)), where, 
neglecting the small spin-orbit coupling, the spin SU(2) 
symmetry is assumed to be present, and so in the sin- 
glet channel M a p = T a <7 v 8 a p while in the triplet channel 
Ma/3 = T u a u a a /3. Similar expressions hold for N\ p and 
/i, v = 0, 1, 2, 3; (to or <ro = H-2). When the quartic terms 
are integrated over r, the symmetries of the lattice re- 
quire that the product of two IR's to which M and N 
belong is the trivial IR. Further simplification comes from 
the fact that all the spin triplet-triplet terms can be writ- 
ten in terms of singlet-singlet terms using Fierz identities 
21, 3l|. Therefore, the most general contact quartic in- 
teraction terms preserving the time reversal symmetry, 
lattice symmetry, and spin SU (2) symmetry are 
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FIG. 1: Schematic zero temperature phase diagram for ABC 
stacked trilayer graphene at the neutrality point, g is the 
strength of the electron-electron interactions and A the con- 
duction and valence band mixing at the symmetry points ±K 
(see Eq. ((2|). In the non-interacting case A c i, C 2 mark the 
quantum phase transitions between distinct, but symmetry 
equivalent, states 3 + , 9~, and 3~. Their schematic contours 
of constant energy near the K point with Dirac points (red 
circles) are shown. The (blue) phase boundary between 3~ 
and 9~ terminates at two tricritical points Pi and P-i- At any 
finite g the direct transition between 3 + and 9~ is avoided 
and replaced by transitions into a symmetry broken phase, 
whose width becomes exponentially small as g — > 0. 



where g lu = g 2v = 9ku, for each v, and ~g x = g 2 = 5k, 
and we inserted the factor of 87t/|to*|, which makes <?'s 
dimcnsionless, for convenience. The total number of in- 
dependent, symmetry allowed, couplings is 15. We dis- 
cuss the specific model choices of g's further in the text. 

To find the phase diagram, we evaluate the partition 
function Z = Tr [e~ ff / T ] using coherent state path inte- 
grals and the Wilson renormalization group (RG) proce- 
dure [32| which, for BLG, has been spelled out in detail 
in Ref. |22J. As the modes within the thin momentum 
shell 1 — d£ < |k|/A < 1 and any Matsubara frequency 
uj n = (2n + 1)ttT are integrated out in powers of small 
<?'s, the remaining modes, V ; kcr( a; ™) ana * their momenta 
are rescaled such that the new |k| cutoff is again A and 

we choose the coefficient of ^—fCrit? in h, lg to remain 

2m* 1 z k 

constant. To 1-loop order we obtain the following RG 
flow equations 



dint 

dS 
df 



dlnr 

_9 — 9 

dt ' 



9^ 

C^oit^S^re) 



(4) 
,(5) 



f = EE A^g*3«* (aP \t t ,St,rt) , (6) 



a/3 kXk,' 
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FIG. 2: T — phase diagrams for the three interaction mod- 
els for TLG studied here, (a) goo-only, (b) weak coupling 
Hubbard model, (c+d) Screened Coulomb interaction. 
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QMfaStM), (7) 



where t e=0 = T/^, 5 e=0 = A/^t and for simplicity 
we set Ri = R 2 = R in Eq. ((2]) giving r e=0 = (2m* A) R. 
The expressions for the coefficients C^, C M , A, A, and 
A's are too long to present here, but they are given in 
the Supplementary Material (SM). $ functions are 



$o (t, f,r) = |£ ^Tr [G^Mta,] 



(8) 



-E 

16 ^ 



^(TrlGA^K)^^]) 2 ^) 



In the above, the one particle Greens function, given by 

-F ht 2u , is evaluated at 



G k 1 (w„) = -icj„ + h^ lg + h 
the circle |k| = A parameterized by tp. 

While the details of the flow described by the above 
equations must be solved numerically, the general trends 
can be readily understood. From Eq. = t£ =0 e 2i 

and r£ = ri = §e~ l . Starting from some fixed initial values 
of g's and sufficiently high temperature t£ =0 , the expo- 
nential decay of the <£>'s with growing £ quickly halts the 
growth of the magnitude of the g's, which saturate to fi- 
nite values as £ — > oo. Such a flow corresponds to a phase 
with no spontaneous SB. Holding the values of the initial 
g's fixed and reducing t^-o, there is an increase in the 
values of £ where <&'s halt the flow, and the limiting mag- 
nitude of the g's increases as £ — > oo. Whether the g's 
eventually diverge upon further lowering of te=o depends 
on the initial values of the remaining parameters. The 
highest ti = o at which g's diverge as £ — > oo is identified 
with the critical temperature t c . At t c , the growth of the 
magnitude of the g's is precisely balanced by the decay in 



the <E>'s. We do not extend our calculations below t c . For 
the parameters corresponding to phases 3 + , 3~ and 9~ 
in Fig.Q] the g's remain finite as £ — > oo even if t£ = o = 0. 

Since <I>o ^ for 5 = 0, a finite S is generated under RG 
even if we set the mixed couplings g^ = as we assume 
later on. This is similar to the 1-loop renormalization of 
the mass (or transition temperature) in bosonic cf) 4 theory 
[33l ] . and is ultimately responsible for the slanted shape 
of the phase boundaries in Fig. [T] near small g. 

To extract physical information about the nature of 
the SB we calculate the ordering susceptibilities, Xi> 
for a large number of excitonic, i.e., particle- hole (p-h), 
and superconducting, i.e., particle-particle (p-p) channels 
C. The C=p-h considered are J2 a p I d 2 ^ a {r)O l a/3 'ip l 3(r), 



where O l — T^a v 8 a p for singlet, and % 



for triplet channels. The C=p-p considered arc 
5 Ea/3 / d 2 Ytp a {Y)6 l al3 i}} {Y) + h.c. The Xi can be cal- 
culated using the methods detailed in Ref. 22 and for 
TLG it is given by 

poo 

/ A£^% llj) e^^\t ll 5 tl n).{lQ) 
Jo 



8tt 



where 



n?(*) = B5^^d€V(€')$^ V ')(^,^,r,0,(ll) 



and constants o^^) and B i '/^ v ,s are given in the SM. 
The above equations hold in the absence of the mixed 
couplings, i.e., = 0. This holds at the bare level in 
any of the model cases we studied in detail. Moreover, 
if absent, g^ is not generated by Eq. (|6]). In each case 
our numerical evaluations of the RG flow equations and 
the Xi ' s benefited from additional analytical study of the 
flow at asymptotically large £ near t c . We were able to 
enumerate all stable flow trajectories (rays) and deter- 
mine the instabilities and exponents for each of these. 

Model 1: Fig. [5^ shows the phase diagram in the for- 
ward scattering limit, i.e., only <?oo(^ = 0) 7^ 0. The 
system exhibits an excitonic (p-h) instability in the spin- 
singlet T3CT2 channel. This SB term belongs to A^ g IR: it 
breaks lattice mirror symmetries and C' 2 rotations (180° 
rotations about the axes perpendicular to principal 3-fold 
axes). It does not break the 3— fold rotational symmetry 
about the principal axes 34 1, inversion, or time rever- 



sal symmetry. We call this mirror-breaking (MB) phase. 
The spectrum inside the MB phase is gapless with DP's 
rotated in the same direction at each ±K valley. Our cal- 
culations indicate that MB is very sensitive to the pres- 
ence of even very small back scattering terms, 5x3/300 
and pko/3oo ^ 4 x 10~ 3 , in which case gapped phases are 
preferred (See Fig. [3]). 

Model 2: The results for the Hubbard model are pre- 
sented in Fig. ^Bp. In this case the initial couplings in Eq. 
© arc 503 = 3oo ~ U, g K o = 9K3 = ^9oo and all other 
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couplings vanish init ially . Due to the SO (A) symmetry of 



the Hubbard model [22|, |35M37| , the phase boundary has 
the same shape for repulsive and attractive on-site inter- 
action U. In the former case, U > 0, we find the lead- 
ing excitonic instabilit y to be Oi = I2C3CT, i.e., a layer 
antiferromagnet (AF) [23[ whose electronic spectrum is 
gapped. AF breaks time reversal and inversion symme- 
try. In the latter case, U < 0, the leading instability is 
found to be equally strong for Oi = I2O3I2 i.e., a sponta- 
neously layer-polarized (LP) state and Oi = tx\i(j% i.e., 
an s ++ -supcrconducting state (SC). The LP breaks in- 
version and the SC the charge U(l) symmetry. Their 
spectrum is also gapped. Just as in the case of BLG, 
a mapping connecting these flows can be found 22 1. At 
the transition into the AF phase, all g's diverge along the 
AF stable ray (see the SM). Due to the SO (4) symmetry, 
the asymptotic flow corresponding to the LP+SC insta- 
bility can be obtained from the AF stable ray. This is an 
unstable fixed ray in the g— space and is a consequence 
of fine tuning. An infinitesimally small interaction that 
breaks the SO (A) symmetry results either in LP only or 
SC only instability. Numerically, we see this as a diver- 
sion of the flow toward one of two distinct stable fixed 
rays corresponding to these phases. 

Model 3: Finally, we consider cases of interaction most 
likely to capture the situation in the experiments where 
e-e interaction is screened due split-off bands and gates: 
backscattering is finite, but small compared to goo- This 
makes c/K3 = .9ko- In the phase diagram, Fig. [2}:, we 
set 3ko/<?oo = 0.1 as a representative value. Changing 
this ratio does not induce qualitative changes as we show 
later in Fig. [3] When g o > the leading instability is 
towards AF, the flow following the same stable ray as in 
the case of the Hubbard model. For goo < 0, the lead- 
ing instability is towards s ++ -SC. The phase diagram in 
the case of attractive back scattering (here represented 
by 5ko/<?oo = —0.1) is shown in Fig. [2ji. For the ma- 
jor part of the transition line, we find that the model 
undergoes a transition to a SB phase with Oi = T^a^a, 
i.e., (gapped in bulk) quantum spin Hall phase (QSH) 
[H]. Only in a small portion of the phase boundary 



mi 



< 0.04) facing the 3 + phase the leading instabil- 
ity changes to Oi = T1CT3I2 or Oi = T2C73IL2, the two 
order parameters being equivalent by lattice symmetry. 
This instability leads to the layer polarization density 
wave (LPDW), a state where both top and bottom lay- 
ers exhibit a charge density wave, but the local density of 
charge has opposite sign on the two layers. The spectrum 
of the LPDW phase is gapless. 

It is easy to understand why spontaneous rotational SB 
phases are unfavorable in TLG. Because the low energy 
degrees of freedom in TLG reside on the sites directly 
above/below each other in the two outer layers, the 120° 
rotations have no effect on the valley-layer matrices T^a v . 
Therefore, any rotational SB order parameter operator 
must involve at least one power of k, making it less rel- 
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FIG. 3: The influence of back scattering on the transition 
temperature and symmetry breaking channel. The color 
of isotherms corresponds to the leading instability: mirror- 
breaking (cyan), anti-ferromagnet (red). In 3* regions no 
symmetry breaking occurs, even at T = 0. Inset: a TLG 
lattice with the AF order. 



evant under RG, and therefore less likely to condense, 
than the k-independent operators considered here. This 
is unlike in BLG, where the two sites are horizontally 
displaced, some t^<j v transform non-trivially under rota- 
tions, and such states can be favored ID, 22|. 

In Fig. [3] we examine the influence of back scattering 
on the transition temperature and the leading instabil- 
ity. We choose to use non-interacting Hamiltonian pa- 
13 that yield A = 0.106/a, r = 4.19, and 
as only for 8 < 0, the TLG phase dia- 



ramcters 
5q = -0.071, 

gram may exhibit reentrant SB phase (see Fig.[2ji). The 
relative strength of back to forward scattering, (?ko/.9oo, 
varies from to 0.1, therefore connecting the left and 
right edges of this plot to Figs. [2^ and c, respectively. 
With no back scattering and at finite T, the critical in- 
teraction goo at which the MB phase sets in must be 
fairly large. A presence of non-zero gKo quickly sup- 
presses this value, implying that the back scattering helps 
to enhance T c . The MB instability is suppressed by back 
scattering, too, with the AF order preferred whenever 
SKo/ffoo ^ 4x 10~ 3 . This result sug gest s that the gapped 
state observed in TLG experiments [24| is most probably 
the AF state (Fig. [3] inset). 

The authors acknowledge discussions with B. Roy in 
the early stages of this work and thank R.E. Throck- 
morton for a critical reading of the manuscript. This 
work was supported by the NSF CAREER award under 
Grant No. DMR-0955561, NSF Cooperative Agreement 
No. DMR-0654118, and the State of Florida. 



[1] B.I. Halperin and T.M. Rice, Solid State Phys. 21, 115 
(1968); Rev. Mod. Phys. 40, 755 (1968). 



5 



[2] K.S. Novoselov, A. K. Geim, S.V. Morozov, D. Jiang, 

Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A. A. Firsov, 

Science 306, 666 (2004). 
[3] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. 

Katsnelson, I.V. Grigorieva, S.V. Dubonos, and A.A. 

Firsov, Nature 438, 197 (2005). 
[4] Y. Zhang, Y.-W. Tan, H.L. Stormer, and P. Kim, Nature 

438, 201 (2005). 
[5] R.T. Weitz, M.T. Allen, B.E. Feldman, J. Martin, and 

A. Yacoby, Science 330, 812 (2010). 
[6] J. Martin, B.E. Feldman, R.T. Weitz, M.T. Allen, and 

A. Yacoby, Phys. Rev. Lett. 105, 256806 (2010). 

[7] A.S. Mayorov, D.C. Elias, M. Mucha-Kruczyriski, R.V. 
Gorbachev, T. Tudorovskiy, A. Zhukov, S.V. Moro- 
zov, M.I. Katsnelson, V.I. Fal'ko, A.K. Geim, and K.S. 
Novoselov, Science 333, 860 (2011). 

[8] J. Velasco Jr., L. Jing, W. Bao, Y. Lee, P. Kratz, V. 
Aji, M. Bockrath, C.N. Lau, CM. Varma, R. Stillwell, 
D. Smirnov, F. Zhang, J. Jung, and A.H. MacDonald, 
Nat. Nanotech. 7, 156 (2012). 

[9] A. Veligura, H.J. van Elferen, N. Tombros, J.C. Maan, 
U. Zeitler, and B.J. van Wees, Phys. Rev. B 85, 155412 
(2012). 

[10] F. Freitag, J. Trbovic, M. Weiss, and C. Schonenberger, 
Phys. Rev. Lett. 108, 076602 (2012). 

[11] W. Bao, J. Velasco Jr., F. Zhang, L. Jing, B. Standley, 
D. Smirnov, M. Bockrath, A. MacDonald, and C.N. Lau, 
Proc. Natl. Acad. Sci. USA 109, 10802 (2012). 

[12] W. Bao, L. Jing, J. Velasco Jr, Y. Lee, G. Liu, D. Tran, 

B. Standley, M. Aykol, S.B. Cronin, D. Smirnov, M. 
Koshino, E. McCann, M. Bockrath, and C.N. Lau, Nat. 
Phys. 7, 948 (2011). 

[13] E. McCann, V.I. Fal'ko, Phys. Rev. Lett. 96, 086805 

(2006), and references therein. 
[14] M. Koshino, and E. McCann, Phys. Rev. B 80, 165409 

(2009). 

[15] K.S. Novoselov, E. McCann, S.V. Morozov, V.I. Falko, 
M.I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, A.K. 
Geim, Nat. Phys. 2, 177 (2006). 

[16] L. Zhang, Y. Zhang, J. Camacho, M. Khodas, and I. 
Zaliznyak, Nat. Phys. 7, 953 (2011). 

[17] In principle, sufficiently large symmetry preserving terms 
may even completely remove the three DF's by moving 
them from K accross the Brillouin zone and annihilating 
them with DF's from — K along the symmetry lines. 

[18] J. Jung, and A.H. MacDonald, arXiv:1208:0116vl (un- 
published) . 

[19] O. Vafek, and K. Yang, Phys. Rev. B 81, 041401 (2010). 
[20] Y. Lemonik, I.L. Aleiner, C. Toke, and V.I. Fal'ko, Phys. 

Rev. B 82, 201408 (2010). 
[21] Y. Lemonik, I.L. Aleiner, and V.I. Fal'ko, Phys. Rev. B 

85, 245451 (2012). 
[22] V. Cvetkovic, R.E. Throckmorton, and O. Vafek Phys. 

Rev. B 86, 075467 (2012) . 
[23] M.M. Scherer, S. Uebelacker, D.D. Scherer, and C. Hon- 

erkamp, arXiv: 1208. 5323 (unpublished). 
[24] W. Bao, L. Jing, J. Velasco Jr, Y. Lee, G. Liu, D. Tran, 

B. Standley, M. Aykol, S.B. Cronin, D. Smirnov, M. 

Koshino, E. McCann, M. Bockrath, and C.N. Lau, Nat. 

Phys. 7, 948 (2011). 
[25] J.C. Slonczewski, and P.R. Weiss, Phys. Rev. 109, 272 

(1958); J.W. McClure, Phys. Rev. 108, 612 (1957); 
[26] B. Partoens, and F.M. Peeters, Phys. Rev. B 74, 075404 



(2006). 

[27] F. Zhang, B. Sahu, H. Min, and A.H. MacDonald, Phys. 
Rev. B 82, 035409 (2010). 

[28] In the absence of particle hole symmetry 9~ would cor- 
respond to a compensated semimental with particle and 
hole Fermi surfaces. 

[29] The order of the transition refers to the Ehrenfest classi- 
fication, i.e., it is labeled by the lowest derivative of the 
free energy with respect to the tuning parameter A which 
is discontinuous at the transition. 

[30] S. Latil, and L. Henrard, Phys. Rev. Lett. 97, 036803 
(2006). 

[31] O. Vafek, Phys. Rev. B 82, 205106 (2010). 

[32] R. Shankar, Rev. Mod. Phys. 66, 129 (1994). 

[33] H. Kleinert and V. Schulte-Frohlinde Critical Properties 

of ^-Theories (World Scientific Publishing, River Edge, 

NJ), (2001). 

[34] Therefore the mirror symmetry breaking state should not 

be confused with Nematic. 
[35] C.N. Yang, Phys. Rev. Lett. 63, 2144 (1989). 
[36] S.C. Zhang, Phys. Rev. Lett. 65, 120 (1990). 
[37] C.N. Yang, S.C. Zhang, Mod. Phys. Lett. B 4, 759 (1990). 



6 



SUPPLEMENTARY MATERIAL 



SYMMETRIES IN ABC TRILAYER GRAPHENE AND IRREDUCIBLE REPRESENTATION 

CLASSIFICATIONS 



The space group of ABC trilayer graphene (TLG) lattice symmetries is P3ml. The lattice is invariant under 
translations by any vector R = m + ib\i + m_i&_i, where m±\ are integers, and b±\ = a\/3/2x ± 3a/2y are the unit 
vectors of the lattice. At T-point, the point group of lattice symmetries is D 3t j; at other high symmetry points of 
relevance to graphene, ±K, the point group of lattice symmetries is D3. The character tables for these two point 
groups are given in Table I. 



D 3d 


e 


2C 3 


3C 2 


i 


2S 6 


3<7<i 


Alg 


1 


1 


1 


1 


1 


1 


Mg 


1 


1 


-1 


1 


1 


-1 


Eg 


2 


-1 





2 


-1 





Aiu 


1 


1 


1 


-1 


-1 


-1 


A 2u 


1 


1 


-1 


-1 


-1 


1 


E u 


2 


-1 





-2 


1 






D 3 


e 


2C 3 


3C' 2 


Ai 


1 


1 


1 


A 2 


1 


1 


-1 


E 


2 


-1 






TABLE I: Character tables for groups D 3c j and D 3 



In addition to obeying the lattice symmetries, we want our Hamiltonian, Eq. (1) in the main text, to be invariant 
under the time reversal, and under the group of global SU(2) spin rotations. The SU(2) spin rotation invariance 
is implemented by taking a general Dirac bilinear 1 i } aMap' l l'p> and noticing that it is a spin-singlet when 

M = T^avdap, while it is a spin-triplet for M = t^ovc?^. The non-interacting part of the Hamiltonian may contain 
only spin-singlet Dirac bilinears, while the interaction part may have either a product of two spin-singlet or a scalar 
product of two spin-triplet bilinears. 

The tensors formed by the remaining two matrices, t^ov, are classified according to the irreducible representation 
(IR) of the space group they transform under: 

A ig,T '■ 3I 4, lo"i; A+ g r : t 3 <t 2 ; A~ g T : r 3 cr 3 ; A~ u T : t 2 1,t 3 cti; A+ u F : ta 3 ; A~ u F : 1<7 2 ; 

A t±K '■ ( T i ± * T 2)1, (n ± ir 2 )ai; A+ ±vi : (ri ± ir 2 )cr 3 ; A 2,±K '■ ( T i =•= iT i)o2- (12) 

The superscripts in the names of IR's denote whether an M component is even (+) or odd (— ) under time reversal. 

The momentum operators and symmetric tensors they form (up to third order) are classified according to the IR's, 
too: 

A+ r :|k|° = l, \k\ 2 = kl + k 2 y , A" r : |k| 3 cos3^ k = Mfc*-3fc.y); A~ u r : |k| 3 sin3^ k = k y (k 2 y - 3*£); 
E+ r : |k| 2 e ±42 ^ = {k 2 x - k 2 y ) ± i{2k x k y ); E~ r : |k|e ± ^- = k x ± iky, |k| 3 e ±l ^ = (fc 2 + k 2 y )(k x ± iky). (13) 

All Dirac bilinears in TLG belong to one-dimensional IR's. Therefore, they are allowed to appcr in a term of the 
non-interacting Hamiltonian, Eq. (2) in the main text, only if multiplied by the momentum tensor belonging to the 
same IR as the bilinear. 



/^-FUNCTION COEFFICIENTS 

Here we present the expressions for coefficients appearing in Eqs. (5-7), (10), and (11) in the main text. 
The coefficients appearing in the flow of <5, Eq. (5) in the main text, are 

C MV = 2Tr ((tritr.) 2 ) - 32^ ^i, and C M = 8 - 32^ - (14) 

Each coefficient A in Eq. (7) in the main text has contributions due to random phase approximation (RPA), vertex 
(V), and ladder, particle- hole (Lph) and particle-particle (Lpp), types of diagrams 

akX.k'X' ■kA,k'A'(RPA) ,kA,k'A'(V) <KA,K'A'(Lph) ■kA,k'A'(Lpp) 
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The same kind of decomposition holds for all A's and A's, appearing in Eqs. (6) and (7) in the main text respectively. 
The coefficients A are therefore the sum of the following 



^jkA,k'A'(RPA) 
^kA,k' A' (V) 



j^k\,k' A' (Lph) 



2<5 MK (5 l/ A^ K "5^A'Tr ( (t^Tq) 2 J Tr ( (<j v ap) 2 



1 



5p K 6„\Tr (T K T K 'T a T K T a T K >) Tr (ctacta'CT/jo-act/jcta') - 

<W<WTr (t k / ^~K^~a T~k,' 7~a7~K ) Tr ((T\XT\CTp(T\>a/3Cr\) 



A 



kX.k' A' (Lpp) 



16 

32 



Tr {TpTnTaTK') Tr {cr v cr\crpcr\> ) Tr (t^ t u t k ) Tr (a v aya/3<Jx) ■ 



(Tr ir^,T K T a T K i ) Tr {a v cr\crp(T\')) + (Tr {T^T K ,T a T K ) Tr {a v a X 'apax)Y 



(16) 

(17) 
(18) 
(19) 



The prefactor P Q( g in the last equation equals +1, for (a/3) = (01), (30), or (32), and it equals —1, for (a/3) = (00), 
(02), or (31). One needs to calculate only coefficients with these (a/3), since is non-zero only when a = or 3, 

and (3 = 0, 1 or 2. This is true for all As, As, A's, a's, and B's, since these coefficients are always summed together 
with $( Q/3 ). 

The coefficients A are given by 



^sr = 2<w<wt>((^) 2 ) 



^oTr ( (0-10-/3) ) +2c5„i 



1«,k'(V) 



^k,k (Lph) 



^jU«Tr {T K T K >T a T K T a T R ,) + S^'Tr (T K <T K T a T K ,T a T K ) 



S^qTy ((0-10/3) 2 ) + 25^1 



- -Tr (r^TaT^ ) Tr (r^r re / r a r K ) 
o 



25^^ ( (crio-^) + (Tr {a-ia v Ofi)Y 



A 



k,k' (Lpp) 



Pal 

16 



(Tr {T^T K T a T K >)) + (Tr (Tfj,T K >T a T K )) 



2S u/3 Tr ( (010/3) I + Tr (0-1 ov 0/3) Tr (01 0730V ) 



Finally, the coefficients A are the sum of the following terms 



TkA,k'(RPA) 



1«A,k'(V) 

M(«/3) 



"IkA^' (Lph) 
M(«/3) 



7kA,k'(Lpp) 
MM) 



= 25 /1K ^ re /Tr ((t m t q ) 



2<5ao + <5aiTi- (V1073) 2 ) 



(20) 
(21) 
(22) 
.(23) 

(24) 



--< 2Jp K Tr (T K T K 'T a T K T a T K i) 



2 ( 5 AO + (5A 1 Tr((a 1 o-, 9 ) 2 ) 



<WTr (T K >T K T a T K >T a T K 



--Tr (t m t k t q t k/ ) Tr (T^T K ,T a T K ) 
o 



— ( Tr ( r A' T « T " T «')) 2 



2 + Tr (o-io-aO/jo-io-^jo-a) 

4<5a/3 + Tr (ctio-ao-,3) Tr (o^oa) 
45a/3 + (Tr(o-io-AO-^)) 2 



(25) 
(26) 



(Tr (r^TaT^Y 



4<5a,3 + (Tr (ctio-^o-a))' 



(27) 



The particle-hole (p-h) instabilities are defined through Dirac bilincars, O l = T K <j\8 a p for singlet and O l = T K <j\<j a p 
for triplet channels. In the case of particle-particle (p-p) instabilities, Dirac bilinears corresponding to singlet channels 
are % = t^gao-^, while those corresponding to triplet channels are % = T K a\a^^ with a = 0, 1, or 3. The coefficients 
in the susceptibilities corresponding to one of these channels are 

r K 

(-1) P < 



^ = 2Tr Ur^f Tr ({a x a v f 



(28) 



Tr I O^r^ut) [O l - [O l j 1 (r M a v l) 
2(-l) p -Tr ((t k t^) 2 ) Tr {<t x ct v <t X (jI) . 



(29) 
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For the p-h channels, the coefficient a depends only on k and A and not on whether the channel is singlet or triplet. 
For the p-p channels, the first line in Eq. (|29[) implies that only channels with antisymmetric O l 's can have a non-zero 
susceptibility. This is already assumed in the second line of the same equation, therefore, it is applicable for singlet p-p 
channel susceptibilities only if t k u\ is symmetric, while for triplet p-p channels it holds only if t k <t\ is antisymmetric. 
Otherwise, cr p_p ) = 0. Note that we used rj = r^, since [i can take values or 3 only. 
Using the same notation, in Eq. (11) in the main text we have 

B i%'»'')' / = 2^ 1K (^A<5f Tr ([t^) 2 ^ Tr ((<7, y ov) 2 ) ~ (t k t,j,Th'T k t^t^,) Tr {a\CF v o v *(j\<y v '<J v ) , (30) 
B &T = -^^Tr (o'WV^f (& - (of) (V^1)(W)) 

= - - ^ Ti (r K T^T^T K T^T^)Tr [axa-^cr^,axau'(ru) ■ (31) 

The symbol Sf takes value 1 if O l belongs to a singlet channel and if it does not. Like before, B^ p ~ p ) vanishes if O l 
is a symmetric matrix. The second line in Eq. (|31[) holds only for antisymmetric p-p channels. 



ASYMPTOTIC BEHAVIOUR OF THE $'S AND THE STABLE RAYS 



We can analyse the RG flows exactly in the t — ► oo limit. This is a consequence of the fact that $ function in Eqs. 
(8) and (9) in the main text, depend only on three parameters, Si, t#, and 7^, one of which will diverge faster than 
the others. Since 7^ flows only by its engineering dimension —1, Eq. (4) in the main text, the dominant diverging 
parameter is either ti or Si. That single parameter determines the asymptotic behaviour of the 4>'s. 

In all the flows that we studied numerically, where < t c < 0.1, we always found Si to diverge faster than ti at the 
phase transition. This is the case we present here. These results are applicable to T = flows, too, since there is only 
one diverging parameter, Si, in that case. We concentrate here on the case with no mixed interaction terms, = 0, 
at the bare level, so none are generated during the entire RG flow. 

Under the assumption that lim^oo ti/Si = 0, we have 



$(00) = _2l 
1*1 



1*1 



(32) 



while $0*") = 0(8)~ 2 for all other $^")'s. Also $ = sgn(<5) + 0(5)~ 1 . These properties allows us to write the flow 
equations in the large t limit, and from there find all the possible stable rays along which the couplings diverge. 
A ray in the coupling space is defined by the following unit vector, 



E E^w, 

\ k=0,3,KA=0 



(33) 



where the summation in the denominator goes only over 12 independent coupling constants, thus (?ka = <?ia or g2\, 
for each A, the two being the same according to the symmetry. A stable diverging ray is given by a unit vector y^i/ 
such that it is a fixed point for the asymptotic flow of z^ v \ 



lim dZ » u 



while at the same time the stability matrix, 







0, 



lim , 



d£ 



(34) 



(35) 



must have no positive eigenvalues. 

For each particular stable ray found, it can be shown that, at the transition point, both Si and g^ v 's diverge as 
e (2+vs)£^ w herc r/s follows from Eq. (5) in the main text, 



ijs = lim ^- y~] g^u{t)C^ (ti, Si, re) . 



(36) 
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One is then required to check if the stable ray obtained is consistent with the initial assumption of Si diverging faster 
than tg. If rjs < for any given stable ray, such a solution to Eq. (|34[) is violating our assumption and should be 
disregarded. 

Using the outlined procedure, we find twelve solutions, the coupling constant stable rays, y^, and present them in 
in Table II. 



No. 


yoo 


Voi 


J/02 


2/03 


2/30 


2/31 


2/32 


2/33 


2/K0 


2/K1 


2/K2 


yK3 


1 





■v/31 














, / 5(6+v^35) 

V 6i 


, / 5(6-%/35) 
V 62 














2 





V31 














. / 5(6 — V3B) 
V 62 


. / 5(6+\/35) 
V 62 














3 





~ V31 


. / 5(6 + v^5) 
V 62 


. / 5(6-V55) 
V 62 


























4 





1 

~ %/3T 


/ 5(6 — V35) 
V 62 


/ 5(6+\/35") 
V 62 


























5 


0.045 


-0.3723 


0.1877 


-0.2825 


0.1394 


-0.2341 


0.6095 


-0.3308 


-0.2341 


0.1394 


-0.2825 


0.1877 


6 


0.045 


-0.3723 


-0.2825 


0.1877 


0.1394 


-0.2341 


-0.3308 


0.6095 


-0.2341 


0.1394 


0.1877 


-0.2825 


7 


0.045 


-0.3723 


0.6095 


-0.3308 


0.1394 


-0.2341 


0.1877 


-0.2825 


0.1394 


-0.2341 


0.1877 


-0.2825 


8 


0.045 


-0.3723 


-0.3308 


0.6095 


0.1394 


-0.2341 


-0.2825 


0.1877 


0.1394 


-0.2341 


-0.2825 


0.1877 
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i 

\/l4 








i 

















3+vT5 
2\/l4 


-3+-/T6 
2\/l4 


10 





1 

VTi 








1 

\7w 

















-3+vf5 
2\/l4 


3+v^5 
2\/l4 


11 


i 

2\/38 


5 

2\/38 


\/15 
2x/38 


x/15 
2\/38 


3 

2i/38 


3 

2v / 38 


2\/38 


\/T5 
2\/38 


3 

2x/38 


3 

2n/38 


Vl5 
2\/38 


2\/38 


12 


1 

2\/38 


5 

2x738 


s/lS 
2\/38 


vTI 

2x738 


3 

2x738 


3 

2\/38 


vT5 
2\/38 


\/l5 
2x738 


3 

2x/38 


3 

2\/38 


%/l5 
2\/38 


%/l5 
2\/38 



TABLE II: Twelve stable rays. 



In Table II, for rays Nos. 5-8 are given by their numerical value, but can be obtained in a closed form in terms 
of zeros of a fourth order polynomial. For example, along the stable ray No. 8, the ratio ui = 2/02/2/03 is a particular 
zero (« —0.5427) of polynomial 

= 7127w 8 + 31470w 7 + 61602w 6 + 91050w 5 + 109022w 4 + 91050cj 3 + 61602w 2 + 31470w + 7127. (37) 

This polynomial is symmetric and can be reduced to a fourth order polynomial by a division by a/ 4 and substitution 
z = lu + 1/uj. Therefore, there is a closed form solution for u which is an algebraic number. The other ratios are then 
related to u> as 

2/oo _ 9 + 14w + 9uj 2 
y 03 8(l + a;) 

2/31 _ VKA _ (1 - Uj) 2 - 

2/03 2/03 8(1 + w) ' 2/03 2/03 16(1 + w) ' 2/03 2/03 16(1 + w) 

where S = (33 + 60a; + 70a; 2 + 60a; 3 + 33a; 4 ) 1 / 2 . In phases Nos. 5-7, the values for the ratios are given by a permutation 
of the No. 8 stable ray ratios. 

The properties of the instabilities associated with each stable ray are listed in Table III. We give the name to 
each phase, list the order parameter of the leading instability, the IR that the order parameter transforms according 
to, whether it is even (+) or odd (— ) under the time reversal symmetry (0), if it is gapped or not, and give rj$. 
The instabilities with order parameter transforming according to an IR at T point have a spatially uniform order 
parameter. The two instabilities with the order parameter transforming according to IR's at ±K represent lattice 
commensurate density waves. 

For the cases of interaction considered in the main text we found flows to go into one of the following phases: MB, 
AF, QSH, LPDW, and s ++ -SC . In the case of the attractive Hubbard interaction, we found the flow to go toward 
a ray obtained from the AF stable ray, Table II, phase No. 8, by keeping y i, y 30 , 2/32; 2/33; an d 2/ki the same and 
reversing the signs of all other y^'s. This ray is fixed (Eq. (|3"4")l holds), but unstable and is a consequence of fine 
tuning of the interaction. 



- 5 2/01 _ -3 + 22a; - 3ar + 35 y 30 _ y^a _ 1 + 0/ 

2/03 16(1 + a;) ' y 3 2/03 2 

2 2/32 2/K2 (1 + 3a;) 2 - 5 y 33 y K3 (3 + a;) 2 - H 



No. 


Name 


Oi 


IR 


(-) 


Gapped 


Vs 


1 


Mirror breaking (MB) 




^2fl,r- sin S let 


+ 


N 


1/3 


2 


Quantum anomalous Hall (QAH) 


r 3 a 3 l 


^2<,,r- sm g let 


- 


Y (in bulk) 


1/3 


3 


Inversion breaking (IB) 


lff2 1 


^2 U ,r- sin g let 


- 


N 


1/3 


4 


Layer polarized (LP) 


lffsl 


^L,r- sin S let 


+ 


Y 


1/3 


5 


Triplet mirror breaking (tMB) 


T 3 a 2 ff 


A^ q r -triplet 


- 


N 


0.4622 


6 


Quantum spin Hall (QSH) 


r 3 a 3 B 


r -triplet 


+ 


Y (in bulk) 


0.4622 


7 


Triplet inversion breaking (tIB) 


1<72<? 


^2u,r- tri P let 


+ 


N 


0.4622 


8 


Layer anti-ferromagnet (AF) 


10-31? 


^2u,r- tri P let 


- 


Y 


0.4622 


9 


Interlayer current density wave (ICDW) 


TKT2I Or T2<T2l 


A^" ±K -singlet 




N 


7/16 


10 


Layer polarization density wave (LPDW) 


T1CT3I Or T2<J 3 1 


Aj ±K -singlet 


+ 


N 


7/16 


11 


s++-superconducting (s++SC) 


ril(72 (p-p) 


^t Q ,r- sin g let 


+ 


Y 


7/16 


12 


Interlayer pairing superconducting (iSC) 


ri a-io-2 (p-p) 


^i"«,r- sm g let 


+ 


N 


7/16 



TABLE III: The list of all the possible phases in TLG with their properties. 



